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ABSTRACT 

We investigate general aspects of molecular line formation under conditions 
which are typical of prestellar cores. Focusing on simple linear molecules, we 
study formation of their rotational lines by radiative transfer simulations. We 
present a thermalization diagram to show the effects of collisions and radiation 
on the level excitation. We construct a detailed scheme (contribution chart) to 
illustrate the formation of emission line profiles. This chart can be used as an 
efficient tool to identify which parts of the cloud contribute to a specific line 
profile. We show how molecular line characteristics for uniform model clouds 
depend on hydrogen density, molecular column density, and kinetic temperature. 
The results are presented in a 2D plane to illustrate cooperative effects of the 
physical factors. We also use a core model with a non-uniform density distribution 
and chemical stratification to study the effects of cloud contraction and rotation 
on spectral line maps. We discuss the main issues that should be taken into 
account when dealing with interpretation and simulation of observed molecular 
lines. 

Subject headings: line: formation — line: profiles — radiative transfer — stars: 
formation — ISM: clouds 



1. Introduction 

Star formation is a fundamental process in the universe. Dense, gravitationally bound, 
and genuinely starless cores, which we call here "prestellar cores", are the earliest visible 
precursors of forming stars. The interplay between gravitation and thermal/magnetic pres- 
sure as well as the conservation of angular momentum are all driving forces behind prestellar 
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core evolution. While the thermal structure is determined by dust properties and molecular 
composition, magnetic support is sensitive to the ionization degree. The interaction and 
relative importance of these pro cesses as well as the role of external forces are sti ll not well 
understood (see, e.g. reviews by lMac Low fc Klessenl 12004 ; iBergin fc Tafallal 120071 ). Thus, it 
is important to study the physical and chemical evolution of these cores in detail to reveal 
the underlying physics. 

The interiors of such prestellar cores are well-shielded from interstellar or stellar radi- 
ation, leading to low internal temperatures. Under such conditions, the main component 
of these cores, H_2, is not easily observable. Therefore, we have to rely on indirect meth- 
ods to determine the physical structure of prestellar cores, e.g., on observations of thermal 
dust emission or emission lines from other molecules, like CO, CS and N 2 H + . Observations 
of spectral lines have an important advantage over continuum observations since they also 
carry information about kinematics of the gas. A disadvantage of molecular tracers can be 
that they are only present under certain conditions and can be frozen-out on dust surfaces. 
Numerous studies, including single-dish and interferometric observations along with their 
theoretical analysis, have been performed over the past years, significantly deepening our 
underst anding of the physical and chemical processes in star-forming regions (s ee, e.g. re- 
views of |Myerjll999|; lEyanj[l999|; IdjjnmcejBCiij^^^ However, 
deriving physical properties from molecular lines is a difficult (inverse) problem because of 
the many factors which can affect the line formation. 

Even if we assume that the studied object is spherically symmetric and uniform, we 
have to specify at least 5 independent parameters describing the formation of molecular 
lines. These are density ^(H 2 ), kinetic temperature T^ in , molecular column density iV(mol), 
radial velocity V r , and micro-turbulent velocity V tUT ^. In more realistic models, additional 
parameters should also be considered, e.g., electron concentration and external radiation 
field, with spatial distributions of all the above parameters. 

An exact analytical treatment of this multi-parameter system is, in most cases, impos- 
sible. Optically thin li ne formation can be described analytically to some extent (see, e.g., 
Rohlfs fc Wilson! |2000| ). but it is very difficult to do the same for optically thick lines. Also, 
it is difficult to use analytic methods to solve the inverse problem, i.e., to restore physical 
distributions of relevant parameters from observed spectra because of the non-local and non- 
linear nature of the radiative transfer problem. This is why numerical line radiative transfer 
(LRT) models are commonly used as an exploratory tool. Several reliable numerical methods 
and numerica l codes have been develope d by various groups for this purpose (see reviews by 
Peraiahl l200ll ; Ivan Zadelhoff et al.l |2002| ) . There are also fast approximate numerical tools 



- 3- 



available for the molecular line radiative transfer analyst (jvan der Tak et aL 2007 ). How- 
ever, input data for line modeling include not only physical conditions but also distributions 
of molecular abundances. All these data can be represented with some analytical prescription 
or extracted from chemical and dynamical models. 



As recently outlined by lTsamis et al.l (120081 ). there are two alternative approaches to the 



diagnostics of protostellar objects. The first approach is a systematic study of the influence 
of various factors on emergent spectra in order to facilitate the analysis of forthcoming ob- 
servations. Such studies (based mainly on approximate LRT methods) have already started 
more than 30 years ago with the analysis of the excitation conditions for various molecules 
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approach is a detailed study of an individual source. To infer its parameters, a number of 
models is constructed and the direct LRT problem is solved for each model. Varying the 
model parameters, it is possible to find the combination which provides the best agreement 
between o bserved and modele d spectra (maps ) or their derived parame ters, like, e.g., velocity 



centroids (jWalker et al.lll994j ) or line ratios (Ivan der Tak et al.l 120071 ). Both trial-and-error 



methods or any sophi sticated numeric al algorithm for isolating the "best" set of free param- 



eters can be utilized ( IKeto et al.l 120041 ) . This second approach has already been successfully 



applied in a number of studies, aimed to extract detailed characteristics of pre-stellar and 



struc ture (e.g., 
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Following the first approach, efforts of various authors are mainly concentrated toward 
evolutionary stages later than the prestellar phase. Most cores investigated in detail are 
closer to the formation of a first hydrostatic core, when the collapse is well-developed and 
can be described by Shu or Larson-Penston solutions. This interest may be partly caused by 



the fact that firs t proofs of infall have been reported for Class objects (jWalker et al.l 11986 



Zhou et al.lll993), in which the collapse has already resulted in the formation of a deeply em- 



bedded source. IZhoul (Il992l ) considered spectral differences between Shu and Larson-Penston 
solutions, assuming constant mole cular abundances, and identified four line properties, sup- 
posedly unique to collapsing cores. iRawlings et al.l (119921 ) combined a dynamical description, 
based on the Shu model, with a detaile d chemical model and ide ntified species with broad 
line wings, indicative of infall motions. IRawlings fc Yated (120011 ) coupled a similar chemi- 
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cal and dynamical model to a more realistic radiation transfer model, implemented as an 
approximate A-iteration code, and studied the sensitivity of emergent central optically thin 
and opti cally thick spec t ra to vario us model paramete r s. Th e assumption of isothermality. 



used b y iRawlings et al.l (119921 ) and iRawlings fc Yated (120011 ) , was relaxed in iTsamis et al. 



( 120081 ) in application to a "B335-like" Class object. The sensitivit y of line profiles in thi s 
model to the intensity of the ambient radiation field was studied in iRedman et al.l (120041 ) . 
The influence of various parameters on Class central spectra, with a parti c ular e mphasis on 
turbulence and rotation, was also studied by IWard-Thompson fc Buckley! (120011 ) under the 
assumption of flat abundance di stribu t ions. The role of high angular resolution in studying 
such clouds was in vestigated by lChoil (120021 ). More details about line modeling can be also 
found in review by lEvand (119991 ) . 



Large amounts of data ha ye also been accu mulated for starless cores, most of which are 
presumably "prestella r" (e.g., ISohn et al.l 120071) . together with the theoretical modeling of 
molecular lines. E.g.. iTafalla et al.l (120041 ) investigated the effects of depletion on the line 
intensity. iLee et al.l (120041 ) studied the evolution of line profiles during core contraction based 



on the models where the chemica 



of the cloud. iPavlvuchenkov et al. 



evolu tion is calculated along with dynamical evolution 
d2007h investigated the combined effects of temperature 
and depletion, and De Vries fc Myers] ( 2005 ) investigated the possibility to use approximate 
analytical models to extract velocity gradients from the line profiles. In principle, all the 
mentioned studies of line formation in protostellar objects are also relevant for prestellar 
cores. However, despite the obvious observational and theoretical progress, the interpretation 
of the data is still far from being straightforward, as infall and rotation velocities in prestellar 
cores are both comparable or even smaller than the sound speed. Difficulties in restoring the 
information about structural, thermal, and kinematic properties are also caused by the lack 
of relevant methods for molecular line analysis. 

In this paper, we systematically study the formation of molecular lines in prestellar 
cores and present different tools which can be used to analyze the formation of lines in 
detail. Using the linear molecules CO and HCO + as examples, we model their emission 
by means of non-LTE LRT simulations. The role of hydrogen density, molecular column 
density, kinetic temperature, infall, and rotation is examined and illustrated in the paper. 
Such an analysis may form the basis of a more sophisticated interpretation of observational 
data and can be useful for those who are going to use LRT simulations in their studies or 
want to get a global view of the factors influencing observed line profiles. 

In Section 2 we describe the radiation transfer model used in the paper, introduce 
parameters to characterize line profiles, and provide some basic considerations on the line 
formation in prestellar cores. In Section 3 we show how molecular line characteristics depend 
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on hydrogen density, molecular column density, and kinetic temperature for the parameter 
sample of uniform model clouds. In Section 4 we consider a non-uniform model cloud and 
study the effects of chemical stratification, contraction and rotation on spectral line maps. In 
Section 5 we discuss additional problems related to the LRT analysis. Section 6 summarizes 
the main conclusions of this paper. 



2. Notes on Molecular Line Transfer 



2.1. Equations of Radiation Transfer 



The goal of the LRT simulation is to obtain level populations and to produce appropriate 
molecular line profiles. For that purp ose, one solves a system of equations, describing the 
radiative transfer (e.g., |Peraiahll200ll ). This system includes the transfer equation itself 



as 

and a set of balance equations for level populations 



X 



^A ul + ^{B ul J. 



ul 



C 



ul) 



l<u 



niAiu + 2j ni(B iu J u i + C, 



lu) 



(2) 



1>U 



Here I u is the intensity of radiation, s is the length along the ray, n u and rii are level popu- 
lations, A u i and B u i are the Einstein coefficients; C u i are coefficients of collisional excitation; 
ul indices specify the transition from the upper level u to the lower level I. Equations ([T|) 
and (j2J) are coupled by the emission and absorption coefficients j u and a u , respectivelly: 
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where <fi u i{v) is the line profile function and Q is the spatial angle. The line profile function 
can be expressed as 
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in the approximation of total redistribution over frequencies and a Maxwellian turbulent 
velocity distribution. Here, v u i is the central frequency of the transition u — > Z; v is the 
regular velocity; n is unit vector associated with dQ; and b is a parameter, which is related 
to the kinetic temperature and the most probable value of the microturbulent velocity 
V t by the expression 

^Zk + V 2 (7) 

^mol 

The intensity I v can be expressed in units of the brightness temperature Tb via the Planck 
equation 

2/iz/ 3 1 

Iv — 2 IE ' ( 8 ) 

C gfcT B _ I 

or in the units of radiative temperature Tr via the Rayleigh- Jeans approximation with sub- 
tracted background radiation 

where J^ g is the intensity of the background radiation with the temperature Tb g . The ra- 
diative temperature Tr is commonly used in radio astronomy and we use it throughout this 
paper. 

In order to calculate the emergent profile, one first calculates level populations and then 
integrates Eq. (pT|) to get a spectrum, which can be later convolved with a telescope beam, if 
necessary. For our LRT simulatio ns we use the 2D non-LTE code URAN(IA), developed by 



Pavlyuchenkov fc Shustovl (12004 ) . This code partly utilizes the scheme originally proposed 



and i mplemented in the publicly available ID code RATRAN (jHogerheijde fc van der Tak 



20001 ). The idea of the method is to solve the system of LRT equations with the Accelerated 



A-Iterations method, where the mean intensity J u is calculated with a Monte-Carlo approach. 

We focus on pure rotational transitions and restrict ourselves to linear molecules with 
the simplest level structure (CO, HCO + ) in order to avoid the complexity caused by blends, 
uncertainties in collisional coefficients, etc. All the molecular data needed for the LRT 
simulation are taken from the Leiden Molecular DatabaJ feri.fi.,, «t all B. 



To describe and to analyze the results of the simulations, we use two important param- 
eters of the transition, namely, the optical depth r v and the excitation temperature T ex . The 
optical depth is defined as 

Sl 

t v — J a u ds, (10) 

so 
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where v is the frequency within the line profile, sq and s\ are near and far edges of the 
cloud along the ray, respectively, and a v is the absorption coefficient. The optical depth 
determines how opaque the medium is and how effectively different regions are coupled by 
radiation. The excitation temperature T ex of the transition is defined as 

ni gi f hu ui \ 

exp _ , (11) 

lb u yu (_ r*' 1 ex ) 

where rii and n u are populations of levels I and u, and gi and g u are their statistical weights. 
The excitation temperature indicates how close the medium (transition) is to local thermo- 
dynamical equilibrium (LTE), when T ex = Tk in . The excitation temperature and the optical 
depth can be used to constrain the LRT problem. In the optically thin case the temperature 
of the emergent radiation T-q is proportional to the optical depth and excitation temperature 
of the transition 

T B « rT ex . (f 2) 

In the optically thick case 

T B « T cx (r = 1), (13) 
i.e., Tb is close to the excitation temperature at r ~ 1. 



2.2. Contribution Chart 



Tafalla et al.l (120061 ) implemented the contribution function (CF) as a convenient tool 
to analyze the relative input of a particular core region to the line profile. They applied it to 
the study of line formation and chemical differentiation in starless cores, using the frequency- 
integrated CF. We extend their analysis, taking into account all the velocity channels, and 
investigate in detail the formation of the entire line profile. The contribution function is 
defined as 

F v (l) = e~ T * {l - e~ Arv } S, (14) 

where t v is the optical depth toward the element at distance I, At v is the optical depth 
of this element, and S is the source function of the element. The index v indicates that a 
corresponding value is related to a velocity offset v. 
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Fig. 1. — Contribution chart for HCO + (3-2) spectrum from uniform cloud. (SPEC) Emission 
line profile toward the center of the cloud. (TAU) Spectrum of the optical depth toward the 
center of the cloud. (TEX) Excitation temperature as a function of position along the ray. 
The position, D, is given in units of cloud radius. The value D = — 1 corresponds to the 
cloud face, while D = 1 corresponds to the cloud rear; (MOL) Distribution of the molecular 
concentration as a function of position. (CF) Contribution function for the given position 
and velocity offset. The projection of regular velocity onto the line of sight is shown with a 
white line. (ICF) Contribution function integrated over the velocity. 
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We find it convenient to combine the contribution function with some other auxiliary 
plots into a common contribution chart, which is shown in Fig. [TJ In this Figure we present 
results of HCO + (3-2) line simulations for a static uniform spherically symmetric cloud having 
T kin = 10 K, ti(H 2 ) = 10 5 cm" 3 , iV(HCO + ) = 10 13 cm" 2 and Kurb = 100 m/s, which are 
typical parameters of prestellar cores. 

The emergent spectrum for the ray passing through a core center is shown in panel SPEC. 
Panel TAU contains the corresponding line optical depth. Distributions of excitation tem- 
perature and molecular number density along the ray are shown in panels TEX and MOL. 
The central part of the chart is occupied by a plot depicting the CF for the ray as a function 
of position and velocity (CF). 

From this chart, one readily sees that in this particular example radiation of the line 
comes almost exclusively from the nearby part of the cloud. Only in wings there is some weak 
contribution from more remote regions. This is emphasized in panel ICF, where the velocity- 
integrated CF along the ray is shown. We will see later that in more realistic situations the 
picture may not be that straightforward. 

2.3. Thermalization and Critical Density 

Before we start modeling the line profiles under the varios conditions, a brief com- 
ment on the conditions in starless cores is in order. There are two primary mechanisms for 
molecules to be excited in molecular clouds. These mechanisms are collisions with other 
species (primarily with H 2 , He, and e~) and radiative excitation. The efficiency of the first 
process depends on gas density and kinetic temperature. The radiative excitation, in turn, 
depends on the intensity of the ambient radiation field. 

To demonstrate the relative importance of collisional and radiative excitation, we con- 
sider a simple model of a uniform cloud with temperature T^ in and molecular hydrogen 
density n(H 2 ), which is illuminated by isotropic blackbody radiation. The temperature of 
background radiation is set to be Tb g = 2.73 K (first case) and Tb g = 50 K (second case), 
while Tkm is varied between 10 K and 100 K. For our analysis, we assume that the cloud is 
extremely optically thin (r <C 1), i.e. we neglect the self-radiation of the cloud to show the 
effect of radiation and collisions more clearly. We discuss the situation for two molecules 
with low and large dipole moments, namely, for CO and HCO + (note, that their lines, as 
a rule, are not optically thin in real prestellar cores). Excitation temperatures of CO (2-1) 
and HCO + (1-0) as functions of hydrogen density and kinetic temperature are presented in 
Fig. [3 
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Fig. 2. — Excitation temperature of CO(2-l) and HCO + (1-0) as a function of hydrogen den- 
sity and kinetic temperature in the optically thin limit. The temperature of the background 
radiation is T^ g = 2.73 K (top row) and 50 K (bottom row). Different colors correspond to 
different excitation temperatures, as indicated on a scale located to the right of each plot. 
The black area in the HCO + (1-0) plots corresponds to negative excitation temperatures, 
i.e., to an inversion in level populations. The solid white line stands for the critical density 
n cr , see Eq. [15] . while the dashed white line corresponds to the thermalization density n t h- 



It can be clearly seen that the map can be divided into three different parts. If the 
hydrogen density is lower than the so-called "critical density" n cr , 



Aui 



(15) 



n, 



cr 



the molecular excitation temperature is determined by the radiation field. In Eq. (I15p . 
A u i is the Einstein coefficient and C w i are collisional rates from the upper level u to other 
low-lying levels i. At densities lower than the critical density, collisions are not effective to 
excite the molecules. In this region (left part of all plots in Fig. |2j) T cx ps Tb g and, thus, 
does not depend on T^ in . Note that n CT is different for different transitions. In our case, 
n CT ps 10 4 cm -3 for CO(2-l), and n cr ps 10 5 cm -3 for HCO + (1-0). Also, n cr is usually 
greater for upper transitions of the same molecule. 

We also introduce the thermalization density nth, as the density, above which collisional 
transitions define level populations, so that levels are "thermalized" , and T ex k, T^ in (right 
part of all plots in Fig. |2J). More precisely, we define Tl\fa clS 8b density at which the relative 
difference between T ex and Tk; n is 5%. 

At intermediate hydrogen densities, n cr < n(H 2 ) < n^, both collisional and radiative 
transitions effectively populate and depopulate levels, so that the excitation temperature is 
equal neither to T kin , nor to T bg . It can even be negative because of specific ratios between 
collisional excitation rates. The important point, which we would like to make here, is that 
these intermediate densities are just typical densities of starless cores. This emphasizes the 
necessity to use adequate methods for interpretation of observed molecular lines. 



In order to perform a qualitative analysis of line profiles, it is necessary to introduce 
some values to characterize them. There are different ways to describe numerically a given 
line profile with a few parameters. In our study we will use the peak intensity T max , the total 
line intensity J (zero moment), the mean velocity V crAr (first moment), the intensity at mean 
velocity T cntr , and the total line width W (second moment). We put a factor (8 In 2) in the 
expression for W in order to make W equal to HPBW for Gaussian profiles. In addition, 
we introduce the relative strength of the self-absorption dip for symmetric profiles, DIP. 
Definitions for these quantities are given in Table [1] and illustrated in Fig. [31 



2.4. 



Characterization of Line Profiles 
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Table 1: Parameters of line profiles 



Parameter Definition 



Peak intensity, T max max Tr(V) 



+00 



Total intensity, J J T R (V)dV 

—00 

Mean velocity, Kntr J' 1 + ]°T K (V)VdV 



—00 



Intensity at mean velocity, T cntr T R (V cntr ) 



-00 



Total line width, W J (8 In 2) J' 1 J (V - V cntr ) 2 T R (V)dV 

y —00 

Self absorption dip, DIP (T max - T cntr )/T max 



Fig. 3. — Various parameters used to characterize line profiles. 
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3. Molecular lines from uniform clouds 

In the previous section we showed that hydrogen density, gas temperature, and back- 
ground radiation are all responsible for line excitation. In the optically thin medium, the 
LRT problem is in some sense local, i.e., the excitation temperature is completely defined by 
local values of rt(H 2 ), Tk in , and by the global radiation field. However, as soon as the total 
column density of emitting molecules is large enough (i.e., lines become optically thick), the 
radiation from the entire cloud plays a role in the molecular excitation. Formation of line 
profiles in the optically thick regime is a complex non-local problem even in the simple case 
of a static uniform spherically symmetric cloud. 

In general, to describe the structure of such a simple cloud one specifies its radius R, 
temperature, Tki n , micro-turbulent velocity, Vt ur b, hydrogen density, n(H 2 ), and the number 
density of a molecule, n(mol). However, the solution of the LRT problem depends on R and 
ri(mol) only via their product R x n(mol) = JV(moljl 

Therefore, spectra of uniform clouds with different radii and different molecular abun- 
dances but with the same molecular column density (other parameters being fixed) should 
be identical. Values of n(H 2 ), T^n, and iV(mol) are sufficient to get a unique solution of the 
LRT problem in a uniform cloud (at some velocity v). From the point of view of balance 
equations, n(H 2 ) and T kin are needed to calculate collision rates, while iV(mol) determines 
the mean radiation field (in the absence of background radiation). Kinetic temperature and 
micro-turbulent velocity affect the line width. 

To show how basic properties of emergent spectra depend on n.(H 2 ), T^ in , and iV(mol), 
we generate a series of models with different parameter sets and perform LRT simulation 
for each of them. Then, for each model we calculate the spectrum toward the center of the 
cloud and extract the main line characteristics (see Tabled]). As three physical parameters 
control the line formation in this setup, we vary two of them, keeping the third parameter 
constant. We present calculated line characteristics as a function of the varied parameters 
to illustrate their cooperative effects, assuming that all these parameters are independent 
in the considered ranges. Molecular abundance X(mol) = n(mol)/n(H 2 ) and temperature 
Tki n are assumed to be uniform over the cloud. For the analysis, we use the HCO + (3-2) line 
which is commonly used in observations to probe dense parts of starless cores. Models in 
this Section are by no means intended to represent real starless cores. We use them only to 
demonstrate the main LRT effects. 



3 It is even possible to use N(mo\)/W as an independent parameter, where W is the velocity dispersion 
(see, e.g., LucasI 1974 ) 
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3.1. Effects of hydrogen density and molecular column density 

In this subsection we consider results of n(H.2) and iV(HCO + ) variations at a fixed tem- 
perature Tki n = 10 K. The cloud is static but has a micro-turbulent velocity Vt ur b=100 m/s, 
which is a typical value for starless cores. The basic characteristics of an emergent spectrum 
of HCO + (3-2) toward the cloud center are presented in Fig. HI 
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Fig. 4. — Greyscale maps of the line parameters for a uniform cloud, a) Peak intensity, T max , 
b) relative depth of the self-absorption dip, DIP, c) total line width, W, and d) logarithm 
of optical depth, lg(r) of the emergent HCO + (3-2) spectrum as a function of n(H2) and 
iV(HCO + ). The kinetic temperature is fixed at 10 K. 
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When both hydrogen density and molecular column density are high (Fig. H^, top right 
corner), T max is close to the LTE value (T R « 4.8 K for HCO+(3-2)). At low n(H 2 ) and/or 
iV(mol) the intensity is also low. However, this plot demonstrates that the common percep- 
tion of a molecule tracing a specific density is somewhat oversimplified. The comparison of 
Fig. @Ji and Fig. Hb shows that the HCO + (3-2) profile can be both bright and self-absorbed 
even if the gas density is much lower than n cr , provided the molecular column density is 
high. 

Of course, this situation may not be realistic f or HCO + (3-2) in p restellar cores as at 



typical HCO + abundances of about few times 10 -9 ( Ohishi et al.lll992l ) the column density 
of 10 14 cm -2 and the hydrogen density of 10 4 cm -3 correspond to a linear extent of about 
3 pc. But, it can be relevant for other molecules and/or transitions, e.g. for CO (3-2). With 
a column density iV(CO) = 10 17 cm -2 (which corresponds to an optical depth of about 100), 
density n(H 2 ) = 10 3 cm -3 (below critical), and molecular abundance X(CO) = 10~ 4 , the 
linear extent of the cloud is about 0.3 pc, which is comparable to typical sizes of prestellar 
cores. 

On top right panel of Fig. HJ the relative strength of the self-absorption dip is shown. 
As long as DIP = 0, the profile is single-peaked, otherwise it is double-peaked. The dip in 
HCO + (3-2) line appears when the molecular column density is larger than 10 13 cm -2 . This 
corresponds to an optical depth «1, and a hydrogen density n(H 2 ) < n CT . If the density is 
greater than the critical density, the line has a single peak even at very high HCO + column 
density. The self-absorption profile forms when a gradient in the excitation temperature is 
present. But at high density the excitation temperature is equal to the kinetic temperature 
over the whole cloud (which is assumed to be uniform). Thus, in such a situation we obtain 
a flat-top profile instead of a self-absorbed profile. 

In the bottom panels of Fig. HI the total line width of the spectra and the optical depth 
are shown as functions of n(H 2 ) and iV(mol). Both parameters are mostly determined by 
the molecular column density and only weakly depend on gas density. 

We consider a broad range of HCO + column densities (five orders of magnitude), with 
lower values being not typical for prestellar cores. Our reason is two- fold. First, we intend 
to show the effects of LRT in general. Second, the presented plots can be relevant for 
the isotopomers of HCO + , e.g. for H 13 CO + and HC 18 + which have similar excitation 
parameters but smaller abundances (defined by both isotopic ratios and isotope-dependent 
chemistry) . 
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3.2. Effects of kinetic temperature and hydrogen density 



In Fig. [5] we show the combined effect of temperature and hydrogen density variations, 
keeping the molecular column density constant. To show the complexity of the LRT effects, 
we explore the case of a moderate optical depth, which corresponds to a representative value 
of iV(HCO + ) = 10 13 cm -2 , as follows from the previous su b sectio n. This value is close to 
that inferred for the CB17 globule in iPavlyuchenkov et al.l (120061 ). To test the sensitivity 
of the considered parameters, we vary the temperature in a range, which is wider than is 
inferred for starless cores, however, all the patterns described below are seen at Tk; n < 20 K. 
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Fig. 5. — Line parameters as functions of ^(H^) and Tk; n at a fixed column density 
iV(HCO + ) = 10 13 cm~ 2 for a uniform cloud, a) Peak intensity, T max , b) relative depth 
of the self-absorption dip, DIP, c) total line width, W, and d) optical depth, r, of the 
emergent HCO + (3-2) spectrum are shown. 
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The increase in temperature as well as the increase in hydrogen density leads to an 
increase of the line intensity. At 10 K, the self-absorption dip only appears at densities 
between 10 4 and 10 6 cm -3 . At densities below 10 4 cm -3 collisional excitation is too slow to 
populate levels of the transition, and the optical depth is also correspondingly low. At larger 
densities (> 10 6 ) levels higher than 3 are populated more efficiently, while lower levels are 
thermalized. Because of this fact, population of lower levels decreases, which results in a low 
optical depth for densities above 10 6 cm -3 . The line width is sensitive to the self-absorption 
dip and thus has its maximum at densities between 10 4 and 10 6 cm -3 . When temperatures 
are higher, the density range, at which the dip appears, gets narrower, and the dip itself 
becomes less pronounced, while the line becomes wider, and its optical depth decreases. 

Even though we consider just one particular molecule and particular transition, the same 
diagrams can be obtained for other molecules and transitions. In essence, the diagrams for 
linear molecules (for the same transitions) look similar. The difference is that features of 
the diagrams are shifted along the corresponding axes reflecting different critical densities 
and different column densities per unit optical depth. Such diagrams for more complicated 
models can be used to illustrate the main non-LTE effects of the line formation process and 
to derive parameters of prestellar cores from the observed spectra. 



4. Molecular line formation in starless cores 

Even in uniform clouds and in the absence of regular velocity field, line shapes can be 
quite different from simple gaussians, having self-absorption dips or flat tops due to thermal- 
ization (see above). Of course, when regular velocity fields and chemical differentiation are 
taken into account, line profiles become even more diverse, which further complicates their 
analysis. In this Section we consider the formation of the HCO + (1-0) line in the presence 
of chemical stratification and regular motions, using a simple "chemo-dynamical" model, 
where the chemical evolution is calculated along with the dynamical evolution of the core. 



4.1. Model of starless core 



We assume that a starless core is spherically symmetric with the density distribution 
given by n(r) = no/(l + (r/ro) p ), where p = 2. 4, rp = 2800 AU, rep = 1.35 x 10 6 cm -3 . Param- 



eters are chosen to represent the L1544 core (ITafalla et al.ll2004l ). To c alculate the chemical 
and ki nematic structure of a starless core, we use the model described in lPavlyuchenkov et al. 
(|2006h . 
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We consider models of a static core, collapsing core, and rotating core. To quantify col- 
lapse and rotation velocities, we assume that the model core forms over a time to = 0.8 Myr as 
a result of contraction from an uniform initial configuration with density n ini = 5 ■ 10 3 cm -3 
and rotation velocity at the core edge Vo = 100 m s _1 . The micro-turbulent velocity is 
100 m s -1 in the static model and 50 m s" 1 in the non-static models. Molecular abundances 



are calculated with the same chemical model as described in iPavlyuchenkov et al.l (120061 ). 
Briefly, it is a time-dependent chemical model which includes gas-phase reactions as well as 
the freezing-out of molecules onto dust grains and their desorption back to the gas-phase. 
Surface reactions are not taken into account. Species stick to dust grains with the proba- 
bility of 0.3 for all components except for H 2 and He, for which zer o sticking is assumed. 
We ca lculate the sticking probability for atomic hydrogen according to iHollenbach &: McKee 
(119791 ) . Three desorption mechanisms are taken into account, which are thermal desorption, 
photodesorption, and cosmic ray induced desorption. It is assumed that a cosmic ray particle 
impulsively heats a dust grain to a peak temperatur e T crp , which is close to 70 K for cold 



0.1/iin silicate grain. Data from iLeger et al. I (11985T) are used to approximate the depen 



dence of T crp on the initial dust te mperature (|Semenov et al.ll2004T) . Gas-phase reactions are 



taken from the UMIST95 ratefile ([Miliar et al.lll997l ). We consider the evolution of species 
containing H, He, C, N, O, Mg, Na, Fe, S, and Si atoms. 

Kinetic temperature is assumed to be 9 K. The core is subject to an external UV field 
with the intensity, G, measured in units of the mean interstellar field (lDrainelll978l ). and to 
cosmic rays, which ionize molecular hydrogen with the rate (. The assumed distance to the 
core is 140 pc as believe d to be relevant f or cor es in Taurus cloud complex (which in fact 
maybe closer; see note in lGoldsmith et al.l (120081 )). Using these parameters we perform LRT 
simulations and produce "ideal" (non-convolved) spectral maps and contribution charts for 
the HCO+(1-0) line. 

To illustrate the importance of the cosmic ray (CR) ionization rate and ultraviolet (UV) 
radiation intensity for the chemical structure of the core in Fig. [6] we present abundance 



profiles for HCO + molecules calculated with G = 0.01 and 1.0 and ( = 10 s and 10 
g- 1 a t t = 0.8 Myr. 
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Fig. 6. — Distribution of the HCO + relative abundances in the prestellar core model for the 
different parameters of UV radiation and CR flux. 
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The relative abundance of HCO + strongly depends on both parameters and varies up 
to a few orders of magnitude. The strong UV radiation destroys molecules in the envelope, 
while enhanced CR flux prevents molecules from sticking to dust grains (via CR-induced 
desorption) and enhances HCO + abundance over the whole cloud. 



4.2. Static core 



We first assume that the core is static, i.e., both infall and rotation velocities are zero. 
The calculated spectral map and the contribution chart for the central spectrum are shown 
in Fig. [7] for the case when G = 0.1 and ( = 10~ 18 s _1 . 
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Fig. 7. — Top: Contribution chart for the central spectrum of the non-uniform static core. 
Bottom: Spectral map of HCO + (1-0) for the same core. 
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The contribution chart demonstrates all the aspects of HCO + (1-0) line formation. In 
the inner part of the core the excitation temperature T ex is nearly equal to the kinetic 
temperature while in the envelope it decreases with radius (TEX). This is mainly because 
of the lower density in the envelope, so that collisions are not effective enough to excite the 
level populations. The molecular density distribution is not uniform but has a maximum at 
R/Ro = 0.2 (MOL). In the inner part the core molecules are depleted onto the dust grains. 
The optical depth in the line center is about 20, i.e., the line is optically thick (TAU). From 
the chart, it can be seen that emission in the line center comes from the outer part of the 
core where the excitation temperature is relatively low. The emission in blue and red peaks 
comes mainly from inside, where the excitation temperature is higher. This results in the 
double peaked profile. At the same time, the very center of the core does not contribute to 
the HCO + (1-0) profile even in line wings because of molecular depletion. 



4-2.1. Uncertainties of the chemical structure 



A non-uniform chemical structure of starless cores is of crucial importance for the in- 
terpretation of observed spectral maps. A particular starless core may look very different 
i n lines of various m olecules just because of the different molecular abundance distributions 
( ITafalla et al.l 120021 ) . In general, current models of starless cores, which include chemical 
and micro-physical processes (e.g., depletion and desorption of molecules on dust grains), 
allow us to understa nd qualitatively distribu t ions of simple mole cules like CO, CS, HCO + , 



and N 2 H+ (see, e.g., iBergin fc Tafallal 120071 : lAikawa et al.l 120081 ) . However, it is still diffi 



cult to use such models for the quantitative analysis, mainly because of two reasons. First, 
rates of many chemical and physical processes, involved in dete rmining the abu n dance gra- 
dients, are still yery u ncertain. Here we refer to the papers by IVasyunin et al.l (120041 ) and 
Wakelam et al.l (120061 ) where the uncertainties of chemical reaction rates are studied. Sec- 
ond, external physical conditions such as intensity of the CR and UV radiation fields, are 
often not well constrained, but they strongly a ffect the core chemical structure (see, e.g., 
Farquhar et al.l 1 19941 ; iPavlyuchenkov et al.l 120061 ) . 



Here we want to illustrate the influence of the G and ( values on the emergent spectra. 
We consider a range of "L1544-like" models, but varying G from 0.01 to 1 and ( from 
10~ 18 s _1 to 10~ 16 s -1 . The resulting difference in abundances affects both the intensity and 
the relative strength of the self-absorption dip (Fig. [S]). 
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Fig. 8. — Results of the variation of the UV and CR radiation fields for a non-uniform model 
of a prestellar core, a) Peak intensity, T max , b) relative strength of the self-absorption dip, 
DIP, c) the logarithm of optical depth, lg (r), and d) the logarithm of contributing density 
lgn(cont), (see text) of the line formation for HCO + (1-0) are shown as functions of G and 

c 
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In the bottom right panel we show the density of the line formation, or "contributing 
density", which we define as the density where the integrated contribution function has a 
maximum (Fig. [7J-ICF). This density approximately indicates which densities are actually 
"traced" by the considered transition. The contributing density depends on the UV and 
CR values and can be even smaller than the critical density (cf. Fig. [2]). In the adopted 
parameter space the optical depth of the HCO + (1-0) transition varies from 10 to 1000. 

From the point of view of the line intensity and its optical depth, ( seems to be more 
important than G. The reason is chemistry. In the considered parameter space, T max and t 
almost do not depend on the UV field intensity, because, when G increases, the maximum 
of the absolute HCO + abundance just moves deeper into core, without changing the column 
density much. On the other hand, an increase in ( affects the entire core, enhancing the 
gas-phase molecule abundance and the column density. 

As can be seen from Fig. HH, the optical depth depends strongly on the column density 
and is almost independent on gas density (i.e., excitation conditions). This is reflected in 
Fig. Efc- Peak intensity and the depth of the self-absorption dip depend both on molecular 
column density and H 2 volume density (Fig. Hlai,b) . Thus, the corresponding plots in Fig.[Hh,b 
have a more "diagonal" appearance. The dip gets deeper as the UV intensity decreases, 
because at lower G the maximum HCO + abundance occurs closer to the core edge, i.e., at 
smaller densities. In other words, when G is low, HCO + (1-0) traces lower densities (Fig.|HJl). 

It must be noted that the description used in this subsection is somewhat simplified in 
the sense that it does not take into account the temperature gradient which may appear due 
to enhanced UV irradiation and/or CR flux. Even though these small temperature variations 
are of minor importance for chemistry, they do affect the excitation conditions and, thus, 
further complicate the situation. 

4.3. Collapsing core 

The next step is to abandon the assumption of zero regular velocity in the core. First, 
we assume that the core is contracting, with a maximum infall velocity of about 50 m/s. The 
corresponding spectral map is shown in Fig. [9] for the case with G = 0.1 and ( = 10~ 18 s _1 . 
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Fig. 9. — Contribution chart for HCO + (1-0) emission from the collapsing core with G = 0.1 
and ( = 1CT 18 s _1 . Because of strong depletion, the line peaks mainly form in the outer 
parts of the core, but there is still some contribution of the core center in the blue peak. 
Projection of regular velocity onto the line of sight is shown with a white line. 



-29- 



In contrast to the static model, line profiles show the classical infall asymmetry, with 
the blue peak dominating over the red peak. The map itself is symmetric relative to the 
center. The combination of factors behind the line profile asymmetry is well seen on the 
chart. The contribution chart for the central spectrum (Fig. ^ has a complex appearance, 
which, in general, follows the projected velocity (white curve in Fig. [9]). There is a common 
perception that optically thick lines trace outer parts of the source. However, as we see from 
the chart in Fig. [9j in general, the line profile includes information about different parts of 
the cloud. The integrated contribution function (Fig. [U]-ICF) shows that both near and far 
parts of the core participate in the formation of the central spectrum. As in the static case 
(Fig. [7]), the emission around the dip (at zero velocity) and in the red peak comes from the 
envelope part closest to the observer. However, this region contributes some emission in the 
blue peak as well. The most of the emission in the blue peak comes from the rear part of the 
model core (a hallmark of the infall asymmetry) and also from its center. The gap between 
the near and rear parts of the contribution function (Fig. [9]-CF) at negative velocities is 
caused by depletion of HCO + in the core center. 



Self-absorption spectra toward centers of some prestellar cores (e.g. iLee et al.l Il999 



20011 ) are dominated by red peaks, as if the corresponding core is expanding. The asymmetry 
of the line profile in a particular core can also be caused by an outflow pow ered by an unseen 
central object (protostar). In addition, spectral maps of some cores (e.g. lLada et~al1l2003h 



indicate that their k i nema tics cannot be described as pure contraction, expansion and/or 



rotation. iKeto et al.l (120061 ) suggested that the velocity field in these prestellar cores can be 
reproduced by a model of an oscillating pressure-bounded, thermally-supported object. 

If we, for the purpose of discussion, consider a model of an expanding core by changing 
the velocity sign, we end up with a spectral map and contribution chart for the central 
spectrum, which are just mirror reflections of the same plots for the model of contracting 
core but with the opposite asymmetry. In general, a detailed dynamical model is needed to 
quantify such motions. 



4.4. Rotating core 

Now, we assume that the core does not collapse, but rotate with the axis being per- 
pendicular to the line of sight. The corresponding spectral map and the contribution chart 
for HCO + (1-0) are presented in Fig. [TUJ The line profiles are asymmetric everywhere in 
the map except for the projection of the rotation axis. To the left of this axis, profiles are 
blue-dominated and blue-shifted along the velocity, while to its right they are red-dominated 
and red-shifted. The map is symmetric with respect to the projection of the rotation axis. 



-30- 



As in the case of the contracting core, the appearance of the contribution function follows 
the projected velocity. 
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Fig. 10. — Contribution chart for HCO + (1-0) emission from the rotating core. Projection 
of the rotation velocity onto the line of sight is shown with a white line. 
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Comparing lines from contracting and rotating models, we see that they produce very 
similar profiles in some off-central positions. Therefore, in order to distinguish between the 
models it is necessary to analyze the complete spectral map, not just the central or single 
off-central spectrum. Moreover, in order to extract the information about kinematics of 
the cloud it is of a great importance to analyze simultaneously spectral maps for different 
transitions (optically thin and optically thick) and different molecules. 



4.5. Collapsing and Rotating core 

Finally, let us assume the more realistic case, i.e. that the core is collapsing and rotating 
at the same time. The corresponding spectral line map and the contribution chart for 
HCO + (1-0) and i = 90° are presented in Fig. [Til The spectra map is no longer symmetric, 
neither with respect to the center (as in the case of a purely contracting core), nor with 
respect to the projection of the rotation axis (as in the case of pure rotation). On the left 
side of the map, the effects of infall and rotation are coherent and profiles are strongly blue- 
dominated. On the right side of the map, infall and rotation tend to produce opposide line 
asymmetries (see Fig. l9T fT0]) and their net result is the formation of nearly symmetric line 
profiles. To the left of the rotation axis, profiles are blue-shifted, while to its right they are 
red-shifted, as in the case of pure rotation. 
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Fig. 11. — Contribution chart for the HCO + (1-0) emission from a collapsing and rotating 
core. The projection of the regular velocity onto the line of sight is shown with a white line. 
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Combination of collapse and rotation allows to reproduce the spe c tral fe atures observed 
in a number of molecular cores. In particular, iPavlyuchenkov et al. used a similar 

model for the CB17 prestellar core and derived the infall and rotation velocity by fitting 
the maps of optically thin and optically thick li nes simulteneously. The presence of bot h 
rotation and infall is also seen in the L1544 core (jOhashi et al.lll999l ; IWilliams et al.ll2006l ). 



(.e.g. 



Walker et al. 


1994; 


Chen et al. 


2007) 



5. Discussion 



Starless cores arguably represent the simplest configuration among all stages of the 
formation of a low-mass star. Despite this simplicity, the interpretation of observed spectra 
is not straightforward, especially, in terms of kinematics, because regular motions in these 
objects have velocities which are comparable to or even less than the sound speed. 

Above, we have given some examples of molecular line formation in static and dy- 
namic cores, demonstrating how line profiles are influenced by gas density and temperature, 
molecular abundance, presence of regular motions, and external parameters which drive the 
chemical evolution. Contribution charts show clearly that any optically thick line, generally 
speaking, does not represent some specific part of the core but rather combines contribu- 
tions from different regions. Incidentally, this is what makes it possible to use optically thick 
lines as a diagnostic for core-wide kinematics. On the other hand, while the contribution 
function in general follows the regular velocity profile, emission at a given velocity is not 
unambiguously related to a region which moves with this velocity, bearing information from 
other regions as well because of finite microturbulent line broadening. 

There are some other potential complications which are not accounted for in the present 
paper. In particular, in the models presented above, the gas temperature is assumed to 
be uniform over the core. In reality, the thermal structure of starless cores can be more 
complex. Following the theoretical considerations, the gas temperature is a result of heating 
by cosmic rays, photoelectric heating, collisional exchange with dust, and cooling by atomic 
and molecular lines. The dust temperature, in turn, is mostly controlled by external radiation 
and expected to rise from inside f~ 5 — 7 K) to outside (~ 15 K) (see e.g. lZucconi et al.ll2001 



Evans et al. 


2001; 


Youne et al. 


2004) 



abundance of key species, like CO, C and C + ( IGoldsmithl l200ll ) . Thus, in order to estimate 



the gas temperature distribution one needs to model th e chemical structure, line radiative 



transfer, and the energy exchange between gas and dust. iKeto fc Casellil (120081 ) showed that 
the thermal structure of starless cores may differ depending on their central density. For 
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cores with high central density (n(H.2) > 10 5 cm 3 ), the temperature rises from inside to 
outside, while for less dense cores the temperature distribution is more uniform, which seems 



to be s upported by observations of several prestellar cores. As shown by lPavlyuchenkov et al. 



(120071 ) over-simplified thermal models may lead to wrong interpretations of line data. Thus, 
the thermal structure of the starless cores should be considered in the LRT simulations. In 
fact, the sensitivity of CO lines to the gas temperature in the envelope makes it possible to 
use these lines to constrain the UV part of the radiation field, while the dust RT modelin g 



allows to estimate the overall level of the interstellar radiation field, (see lEvans et al.ll2005l ). 



From the observational side, there is a convolution problem which is related to the finite 
resolution of radio-telescopes. If the telescope beam is comparable to the angular size of an 
object then different regions contribute to the emergent line profile which makes derivation 
of source parameters even more difficult and often ambiguous. 

There are also some uncertainties which are related to the LRT problem itself. One of 
them is the commonly used assumption of the full frequency redistribution (FFR) over the 
line profile which dramatically simplifies equations of the LRT (e.g., the source function does 
not depend on frequency). However, this assumption may not be valid in low density regions 
where collisions are not frequent enough to redistribute the absorbed energy. As a result, 
the line formation may proceed in a way between two extreme modes, namely, between FFR 
and coherent scattering. In the latter case, lines would not have any self-absorption dips 
(but shapes can still be complex because of the regular velocity). 

Another commonly used approach is to simplify the non-regular (turbulent) velocity 
field in cores by introducing the micro-turbulence velocity (Eq. [7]). Such representation 
of the real (unknown) velocity field by the combination of the regular and microturbulent 
velocities also may lead to misinterpretations. For instance, line asymmetry can be consid- 
ered as an evidence of infall (outflow), being actually a reflection of a complex kinematic 
structure due to turbulence or oscillations. Moreover, even if we are not interested in line 
shap es, the microturbu l ent ap proach can fail in reproducing line intensities. As shown by, 
e.g.. iHegmann &: Kegell (120001 ). the turbulent structure of prestellar cores can be rather rep- 
resented by meso-turbulence which again requires the modification of the radiative transfer 
equation. 



6. Conclusions 



In this paper we analyze molecular line formation at conditions typical of starless cores. 
In particular, we show the effect of density, molecular column density, and temperature 
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on line parameters for a sample of uniform clouds. We also consider non-uniform models 
and show the effect of the chemical differentiation, collapse and rotation on the molecular 
spectral maps. We present a chart of line formation which may serve as valuable tool for 
understanding results of LRT simulations. This chart clearly demonstrates which parts of 
the model core contribute to the line profile at each velocity. 

We would like to make the following conclusions: 

• Densities in starless cores fall into the range where level populations are neither radia- 
tively, nor collisionally dominated. 

• Large column densities do not necessarily lead to the appearance of self-absorption 
dips. 

• When the column density is fixed, a specific line can be optically thick only in a range 
of densities, being optically thin at densities both below and above this range. 

• The density which is "traced" by some transition depends on external factors (UV 
field and CR ionization), which shape the molecular distribution. In particular, the 
"traced" density can be lower than the critical density. 

• Rotation and infall may produce very similar spectra and in general can only be dis- 
tinguished by spectral mapping. 

We are grateful to Ted Bergin, Michiel Hogerheijde, Kees Dullemond and Juergen 
Steinacker for useful discussions. We also thank the referee Neal J. Evans II, for valu- 
able suggestions and comments. This research has made use of NASA's Astrophysics Data 
System. DW, and BS are supported by the RFBR grant 07-02-01031. 
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